namespace NBody.Analysis
{
    using System;
    using DataFile;
    
    public class EnergyBinDOS : EnergyBin
    {
        public override void Analyse(NBodySystem s)
        {
            Console.Error.WriteLine("Calculating density of states");
            Console.Error.WriteLine("  Binning in energy with normal spacing");
            // calculate total energies.
            double e_tot;
            Energy.CalculateTotalEnergy(ref s, out e_tot);
             
            // find max, min e
            double emin = ((EnergyParticle)s[0]).Energy;
            double emax = emin;
            foreach (EnergyParticle p in s)
            {
                if (p.Energy > emax)
                    emax = p.Energy;
                if (p.Energy < emin)
                    emin = p.Energy;
            }
            emax += 1e-10;
            Console.Error.WriteLine("  Emin = {0:f6}\n  Emax = {1:f6}", emin, emax);
            
            double delta_e = (emax - emin) / (double)n;
            data = new Table(n, 2, "Density of states (radial orbits)");
            for (int i = 0; i < n; i++)
            {
                data[i, 0] = emin + 0.5 * delta_e + i * delta_e;
                data[i, 1] = 0.0;
            }
            
            double[] r = new double[s.NumParts];
            for (int i = 0; i < s.NumParts; i++)
                r[i] = s[i].Radius;
            
            double[] g1 = new double[s.NumParts];
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < s.NumParts; j++)
                {
                    g1[j] = data[i, 0] - ((EnergyParticle)s[j]).Potential;
                    if (g1[j] < 0.0)
                        g1[j] = 0.0;
                    else
                    {
                        if (orbit == Orbit.Radial)
                            g1[j] = 1.0 / Math.Sqrt(g1[j]);
                        else
                            g1[j] = Math.Sqrt(g1[j]) * r[j] * r[j];
                    }
                }
                if (orbit == Orbit.Radial)
                {
                    // I don't know which is right here ...
                    //data[i, 1] = Math.Sqrt(2.0) * Integration.Trapezoidal(r, g1, 0, s.NumParts - 1);
                    data[i, 1] = 8.0 * Math.PI * Math.PI / Math.Sqrt(2.0) * Integration.Trapezoidal(r, g1, 0, s.NumParts - 1);
                }
                else
                    data[i, 1] = 16.0 * Math.PI * Math.PI * Math.Sqrt(2.0) * Integration.Trapezoidal(r, g1, 0, s.NumParts - 1);
            }
        }
    }
}
            
